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Abstract 

We present a mesoscopic model of the fluid-wall interactions for flows in microchannel 
geometries. We define a suitable implementation of the boundary conditions for a discrete 
version of the Boltzmann equations describing a wall-bounded single phase fluid. We dis- 
tinguish different slippage properties on the surface by introducing a slip function, defining 
the local degree of slip for mesoscopic molecules at the boundaries. The slip function plays 
the role of a renormalizing factor which incorporates, with some degree of arbitrariness, the 
microscopic effects on the mesoscopic description. We discuss the mesoscopic slip properties 
in terms of slip length, slip velocity, pressure drop reduction (drag reduction), and mass flow 
rate in microchannels as a function of the degree of slippage and of its spatial distribution and 
localization, the latter parameter mimicking the degree of roughness of the ultra-hydrophobic 
material in real experiments. We also discuss the increment of the slip length in the transi- 
tion regime, i.e. at 0(1) Knudsen numbers. Finally, we compare our results with Molecular 
Dynamics investigations of the dependency of the slip length on the mean channel pressure 
and local slip properties (Cottin-Bizonne et al. 2004) and with the experimental dependency 
of the pressure drop reduction on the percentage of hydrophobic material deposited on the 
surface - Ou et al. (2004). 



1 Introduction 

The physics of molecular interactions at fluid-solid interfaces is a very active research area with 
a significant impact on many emergi ng applications in material s c ience, chemistry, micro / nano- 
engine e ring, biology and medicine, see lWhitesides fc Stroockl ll200lh . lLion et al\ ll2003lblGad-el'Hakl 
l|l999|) . lHo fc Taifl|l998l) . As for most problems connected with surface effects, fluid-solid inter- 
actions become particularly important for micro- and nano-devices, whose physical behaviour is 
largely affected by high surface/ volume ratios. Recently, due to an ever-increasing interest in mi- 
crofluidics and MEMS (micro-electromechanical system)-based devices, experimental capabilities 
to test and analyze such systems have undergone remarkable progress. 

In this paper, we shall focus on flows in micro-channels, a subject which has recently become 
accessible to sys t ematic experimental studies tha nks to the developments of silicon technology (see 
iTabelinel l)2003|) , iKarniadakis fc Benskokl l|2002l) and references therein) . 

Classical hydrodynamics postulates that a fluid flowing over a soli d wall stick s to th e bound- 
aries, i.e. the fluid molecules share the same velocity of the surface, iBatchelorl l|l989l) : IMassevI 
l|l989|) . This law, and its consequences, are well verified at a macroscopic level, where the char- 
acteristic scales of the flow are much larger than the molecular sizes . The situation c hange s 



drastically at a micros c opic level. Many experiments (iMaure^e^a/J u et al 
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iTrethewav fe Meinhard (|20Q2^I : iBonnacurso et all l|2002l l2003tl: IChoi et all (120031): IZhang et all 
(2002)) and numerical simulations us ing Molecul ar Dynamics llBarraT^TFjocaued lll999alblklThompson fc Troianl 
|l997[) : En"ompson fe Robbind |l99"o(klBocauet fe Barratll)l993(l : ICieplak et all ()200l|) : lThompson fe Robbind 
l|l989() : IPrieziev et all l)2004|V ) have shown evidences that the solid-fluid interactions are strongly 
affected by the chemico-physical properties and by the roughness of the surface. For example, 
water flowing over a hydrophilic, hydrophobic or super-hydrophobic surface, may develop quite 
different flow profiles in micro-structures. One of the most spectacular effect is the appearance 
of an effective slip velocity, V s , at the boundary, which, in turn, may imply a reduction of the 
kinetic energy dissip ation with significant e nhancement of the overa ll throughput (at a give n 
pressure drop), see IWatanabe et all (If 9991); IVinogradova fe Yabukovl (120031): IPit et qZ.I (|2000[) ; 
Baudrv et all (l200lh:ICraig et all ll200lfklZhu fc Granickl |200lll2002|):ICheng fc Giordano! ll2002l) : 



Trethewav fc Meinhartl (|2002^ l ; IBonnacurso et all (|2002j) ; IChoi et aJpOO^ . From the slip veloc- 
ity, one defines a slip length, L s , as the distance from the wall where the linearly extrapolated 
velocity profile vanishes. The experimental and theoretical picture is still under active develop- 
ment. No clear systematic trend of the slip effect as a function of the chemico-physical components 
has been found to date. Slip lengths varying from hundreds of ran up to tens of /jm have been 
reported in the literature. Moreover, controversial claims about the importance of the roughness 
of the surface and of the combined degree of roughness-hydrophobicity have been presented. In 
simple flows roughness is expected to increase the energy exchange with the boundaries, inducing 
a corresponding decrease in the slippage. However, both increase and dec rease of the slip l ength 
as a f unction of the surface roughnes s have been claimed in the literature. IZhu fc Granickl (|2001, 
20021) : IBonnacurso et all l)2002Ll200^1 . From a purely molecular point of view, a critical parameter 
governing the solid- liquid interface is the contact angle (wetting angle). Clean glass is highly hy- 
drophilic, with an angle with water close to 9 = 0° (perfect wetting). Recently ultra- hydrophobic 
surfaces have been obtained which prove capable of sustaining a contact angle with water as high 
as 9 = 177° a value at w hich water droplets are almost spherical on the surface. [Chen et'all (^999); 
iFadeev fc Ca.rthvl fmflj. 

Some authors proposed that the increase in the slippage might be due to a rarefaction of the 
flow close to the wall, a deple ted water r e gion o r vapor layer s h ould e xists near a hydr o phobi c 
surface in contac t with water, iLum et all ()l999|) : ISakurai et all (11998^) : ISchwendel et all l)2003j) ; 
iTvrrell fc Attardl (200^). Recent Molecular Dynamics simulations have also presented some evi- 
dences of a dewetting transition, leading to a strong increa se of the slip length, below some ca pillar- 
ity pressure in microchannels with heterogeneous surface. ICottin-Bizonne et tll ll200fl l20041. The 
physics looks very fragile, depending as it seems on many complicated chemical and geometrical 
details. 

From the numerical point of view, Molecular Dynam i cs (MP) is t he standard tool to s ystem - 
atically investigate the problem, iFrenkel fc sTmH l|l995l) : iRapaportl (|l995j) : iBoon fc Yir3 |l99l|) . 
In MD the solid-liquid and the liquid-liquid interactions are introduced by using Lennard-Jones 
type potential (with interaction energies and molecular diameters adjusted from experiments). 
By changing the interaction energies one may tune the surface tension and, consequently, the 
contact angles. MD also offers the possibility to model the boundary geometries and roughness 
with a high degree of fidelity. The main limitation, however, is the modest range of space and 
especially ti me-scales, which c a n be simulated at a reason able computational time, typically a few 
nanoseconds, iRaoaporH ( 19951) ; iKoolik fc Banavail l|l99ll) . 

The coupling between MD and hydrodynamic modes involves a huge gap of space and time 
scales. Recently, an interesting attempt to reproduce MD simulations of heterogeneous mi- 
crochannels with a continuum mechanical description based o n Navier-Stokes equat i ons and suit - 
able hydrodynamic b ou ndary conditions has been p roposed, (jCottin-Bizonne et 
IPrieziev et all l)2004^ ). ICottin-Bizonne et all l)2004l) show that some of the results obtained by 
MD simulations of a microchannel with a grooved surface can be qualitatively reproduced using 
a Stokes equation for the incompressible flow, in combination with an heterogeneous boundary 
condition, linking the slip velocity parallel to the boundary U||(r) to the stress in the normal 
direction, h: 
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U||(r)=&(r)3 n u||(r) (1) 

where 6(r) is a position dependent normalized slip length mimicking the heterogeneity of the 
microscopic level. The qualitative agreement with the results of MD simulations can be obtained by 
properly tuning the b(r) values. In particular, they show that the dewetting transition observed in 
MD simulations, for some values in the Pressure- Volume diagram, is equivalent to the assumption, 
at the hydrodynamic level, that the boundary surface is made up of alternating strips of free-shear 
(high slip length b(r)) and wetting material (low slip length). 

In this paper, we main aim at filling the gap between the microscopic description typical of MD, 
and the macroscopic level of the Navier-Stokes equations by using a mesoscopic model based on 
the Boltzmann Equation. In particular, we will use a discrete model known as Lattice Boltzmann 
Equation (LBE) with heterogene ous boundary co nditions. 

The boundary condition Maxwell! I|1879|) . arises naturally in a power expansion of the 
Boltzmann equation in terms of the Knudsen number, Kn = A/L, that is, the ratio between the 
mean free path, A, and a typical length of the channel, L. At first order in Kn, one obtains th e 
Navier-Stokes equation w ith the Maxwell boundary conditions above. ICercignani fc Danerill|l963h : 
lHadiiconstatinoul l)2003h . Still, recent experimental results raised s ome doubts on the va lidity of 
this construction above some critical value of the Knudsen number. iMaurer et all (|2004fl . There, 
the authors report that above Kn ~ 0.3 ±0.1 both helium and nitrogen exhibit a non linear 
dependence of the flow rate on Kn which cannot be explained by solving the Stokes equation with 
the first order slip boundary condition Q. For those values of Kn, the flow is in the so-called 
transition regime and it has been shown that the coupling between hydrodynamic equations with 
a second order boundary condition 

U||(r)=6 1 (r)a„U||(r)+6 2 (r)^u || (r) (2) 

is more appropriate to fit the experimental data lMaurer et al\ l|2004l . 

The purpose of our investigation is twofold. First, we aim at developing a model which allows 
a coarse-grained treatment of local effects close to the flow-surface region, without delving into the 
detailed molecule-molecule description typical of MD. Second, we wish to design a tool capable of 
describing fluid motion also beyond the linear Knudsen regime. 

The underlying hope behind the present hydro-kinetic approach, is that the main features of 
the fluid-surface interactions can be rearranged into a suitable set of renormalized LBE boundary 
conditions. This implies that all details of the contact angle, the solid- fluid interaction length, the 
local microscopic degree of roughness, can, to some extent, be included within the local definition 



of effective accommodation fac tors g o verning the statistical interaction s between the mesoscopic 
molecules and the solid walls, ISuccil l|2002() : ISbragaglia fc Succil l)2004|) . In a more microscopic 
vein, one may also describe the interactions betwe en solid-liquid and liquid-liquid populat i ons us 



vein, one may also describe the interactions between solid-liquid and liquid-liquid populations us- 
ing, a mean-field multi- p hase LBE desc r iption , see IShan fc Chenl l)l993l Il994j) : ISwift et al\ ()l995|) : 
IVerberg fc Laddl fcOOOl) : IVerberg et al\ fcOO^ : iKwokl ll20Q3h. Res ults based on these more sophis- 
ticated schemes will be reported in a forthcoming paper. iBenzi et all I^OOfil) . 

The paper is organized as follows. In section JSJl we briefly remind the main ideas behind 
the lattice versions of the Boltzmann equations and we present a natural way to implement non- 
homogeneous slip and no-slip boundary conditions in the model. In section Q we discuss the 
hydrodynamic limit of the LBE previously introduced with particular emphasis on the form of 
the hydrodynamic boundary conditions in presence of slippage. In section Q we present the 
numerical results at various Knudsen and Reynolds numbers, as well as a function of the degree 
of slippage and localization. Whenever directly applicable we compare the results obtained within 
our mesoscopic approach wit h (i) exact results in the limit of small K nudsen numbers obtained in 
the hydrodynamic formalism. IPhilitl lll972albl):lLauga fc Stonel J2003I) (ii) results obtained with a 
microscopic approach using MD simulations, ICottin-Bizonne et al l ( | 2004f) and (hi) recent exper- 
imental results of microchannels with ultrahvdrophobic surfaces. lOu et al\ ((2004^ . Conclusions 
and perspectives follow in section (jSJ). Technical details are given in the appendices. 
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2 Lattice Kinetic formulation 



The Boltzmami Equation describes the space-time evolution of the probability density /(r,v,i) 
of finding a particle at position r with velocity v at a given time t. This evolution is governed by 
the competition between free-particle motion and molecular collisions which promote relaxation 
towards a non- homogeneous equilibrium, whose distribution f eq (p, u), is the Maxwellian consistent 
with the local density, p(r), and coarse grained velocity, u(r). The hydrodynamic variables are 
obtained as low-order moments of the velocity distributions. Infact, the hydrodynamic density 
and velocity are p(r,t) = J dvf(r, v), and u(r, t) — f dvv/(r, v), respectively. The Navier- 
Stokes equations for the hydrodynamic field s are r e cover ed in the limit of small-Knudsen numbers 
using the Chapman-Enskog expansion ICercienanil l)l99l|) . The Boltzmann equation lives in a six- 
dimensional phase-space and consequently its numerical solution is extremely demanding, and 
typi cally hand led by stochastic methods, primarily Direct Simulation Monte Carlo (for a review 
see lBirdl l)l998|) ). However, in the last fifteen years, a very appealing alternative (for hydrodynamic 
purposes) has emerged in the form of lattice versions of the Boltzmann equations in which the 
velocity phase space is discretized in a minimal form, through a handful of properly chosen discrete 
speeds (of order ten in two dimensions and twenty in three dimensions -see appendix A for details). 

This leads to the Lattice Boltzmann Equations (LBE) for the probability density, /;(r, t), where 
r runs over the discrete lattice, and the subscri pt I = 0,N —1 labels the TV discr e te velocities 
values allowed by the scheme, v £ |c n , • • ■ c N _i j.lSuccil f200lh : IWolf-Gladrowl l)2000|k LBenzi et all 
i|l992j) : lc"hen fc Pooler] l|l998|) : lMcNamara fe Zanettil l)l998|) . It is interesting to remark that it is 
sufficient to retain a limited numbers of discretized velocities at each site to recover the Navier- 
Stokes equations in the hydrodynamic limit. In two dimensions the nine-speed 2DQ9 model 
(N = 9) is in fact one of the most used 2d-LBE scheme, due to its enhanced stability iKarlin et all 
(1999). All three-dimensional simulations described in this paper are based on the the 3-DQ19 
scheme (N = 19) (see Fig. ^|in appendix A for a graphical description of LBE velocities in 2d 
and 3d). For the sake of concreteness, we shall refer to the two-dimensional nine-speed 2DQ9 
model, although the proposed analysis can be extended in full generality to any other discrete- 
speed model living on a regular lattic e. We begin by consider ing the Lattice Boltzmann Equation 
in the following BGK approximation. iBhatnagar et all l) 19541) : 

/,(r + c h t + 1) - /,(r,f) = -i (/,(r, t) - f[ eq) (p, u)) + F, (3) 

where we have assumed lattice units 8 X — S t — I. In ©, r is the relaxation time to the local 
equilibrium, which is proportional to the Knudsen number. The explicit expression of the speed 
vectors, c/, of the lattice equilibrium distribution, fi £q \p, u) and of the forcing term F; needed to 
reproduce a constant pressure drop, are described in the appendix A. The hydrodynamic fields in 
the lattice version are expressed by: 

p(r) =£/,(!•); p(r)u(r) =^cj/j(r). (4) 

i i 

Boundary conditions for Lattice Boltzm ann simulations of m i croscopic f l ows have made the objec t 
of much investigation in r ecent years. iToschi fc Succil ll2f)05MNiu et all ll20f)4h : iLim et all (12002^1: 
lAnsumali fc Karlinl ((2002) . In particular, we are interested in studying the evolution of the LBE 
in a microchannel with heterogenous boundary conditions (H-LBE) -the simplest case being a 
sequence of two alternating strip with different slip properties, as depicted in Fig. (JJJ. A general 
way of imposing the boundary conditions in the LBE reads as follows: 

f- k {v w ,t+l) = Y J B- k . l {v w )f l {v w ,t) (5) 
I 

where the matrix Bj. j is the discrete analogue of the boundary scattering kernel expressing the 
fluid- wall interactions. Here and in the following, we use the notation to indicate the generic 
spatial coordinate over the surface of the wall and the indices I, k label the subset of incoming and 



4 



H 




Figure 1: Typical geometry of the microchannel configuration. We have periodic boundary con- 
ditions along the stream-wise, x, and span-wise y directions. The two rigid walls at z = 0, L z are 
covered by two strips of width H and L — H, where L = L x for transversal strips (left panel) 
and L — L y for longitudinal strips (right panel). The two strips have different slippage properties 
identified by the values So and s%. The ratio £ = H/L identifies the fraction of hydrophobic 
material deposited on the surface. Typical sizes used in the LBE simulations are L x — L y = 64 
grid points and L z — 84 grid points. This would correspond, for example, for a fluid like water at 
Kn = 10~ 3 , to a microchannel of height of the order of 100 fim. 



outgoing velocities respectively. To guarantee conservation of mass and normal momentum, the 
following sum-rule applies: 

E%( r ») = L (6) 

k 

Upon the assumption of fluid stationarity, we can drop the dependence on t and write: 

h = Y, B hii^)fi- ( 7 ) 

i 

The simplest, non trivial, application involves a slip function, s(r w ), representing the proba- 
bility for a particle to slip forward, (conversely, 1 — s(r w ) will correspond to the probability for 
the particle to be bounced back). If we focus, for example, on the north- wall boundary condition 
(see Fig. 1)100 in appendix A), the boundary kernel upon the assumption of preserved density and 
zero normal component of the velocity field © takes the form : 

1 - s(r w ) s(t w ) \ ( h\ 

1 / a . (8) 

s(r w ) 1 - s(r w ) J \ U J 

In this language, the usual no-slip boundary conditions are recovered in the limit s(r w ) — > 
everywhere (incoming velocities are equal and opposite to the outgoing velocities), while the 
perfect free-shear profile is obtained with s(r w ) = 1. The formalism is sufficiently flexible to allow 
the study of both spatial inhomogeneity of a given hydrophobic material and/or the effects of 
different degrees of hydrophobicity at different spatial locations. 

The above LBE sche me has been already succ essfully tested in the case of a homogeneous 
slippage s(r w ) = s , Vr,„. ISbragagha &: Succil (|2004^ . In that case, it has been shown (see appendix 
B) that the LBE scheme converges to an hydrodynamic limit with the slip boundary condition 

U||(r w ) = A Kn Idnu^)] + B Kn 3 \%u\\(t w )\, (9) 
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where the parameters A, B can be tuned by changing the degree of slippage, sq and the external 
forcing. In this case, the LBE reproduces the analytical prediction for the slip length, obtained 
by assuming the existence of a Poiseuille velocity profile and, with a suitable choice of A, B in (JSJ, 
one can show that the model is also able to fit th e experimental non-linear dependencies on the 
Knudsen number observed in lMaurer et al\ ()2004(l for nitrogen and helium. 

3 Hydrodynamic limit 

To begin with, we wish to analyze the hydrodynamic limit, Kn — > 0, of the previous LBE models 
with non-homogeneous boundary conditions, as dictated by the space-dependent profile of the 
slip function, s(r w ), at the walls. For the sake of simplicity, we shall confine our attention to the 
continuum limit of zero lattice spacing and time increments, 8 X = 5t — ► 0, c = ^ — ► 1. Starting 
from the discretized equations © one gets for the continuum limit of the LBE: 

dtfi + (c, • V)/, = ~ (ft - fl eq) ) + F t . (10) 

In the following, we shall be interested in the case of stationary, time-independent, solutions 
(small- Reynolds regime). To this purpose, we may formally write the solution of (|1U|) by using the 
time-independent Green's function: 

oo 

/<(r) = £(-ir(r(c ; V)) 

n=0 

Let us notice that defining r — g " £ * (c s being the sound speed velocity), the above expression 
can be interpreted as a formal solution in powers of the Knudsen number. By recalling the 
expression of the hydrodynamic fields Q , it is readily checked that the boundary velocity can be 
expressed as a function of the velocity stress, diUj, at the boundary itself. For sake of simplicity, 
we report here only the first order term (in the Knudsen and Mach numbers) of the expansion 
(see Appendix B): 

U||(r w ) =Kn (-) ^ |9nU||(r w )| + 0(Kn 2 ) (12) 

which is a direct generalization of the result obtained for the case of homogeneous boundary 
conditions ©• The main difference is that, due to the spatial dependence of the stress tensor 
along the wall, subtle non-linear effects may be triggered by the spatial correlation between the 
slip function s(r w ) and the stress at the wall. 

The hydrodynamic equations of motion in the stationary case read as follows: 

r ( u -V)u = -2f + iV.(fpVu) 
I V • {pa) = 

1 u n (T w ) = Kn^) T ^\d n u n (r w )\ + 0(Kn*) (13) 
[ u ± (r w ) = 

where VP contains both the imposed mean pressure drop, F, and the fluid pressure fluctuations. 
In the limit of small Mach numbers <C 1) we may take a constant density p = 1. Let us notice 
that in this limit, the incompressibility constraint V • u = imposes that any non- homogeneity 
of U|| along the wall-parallel direction must be compensated by an equal and opposite gradient of 
the normal velocity u^. This implies that the local velocity profile cannot be of Poiseuille type 
everywhere (uj_ = 0). 

In order to assess the effects of the slip on the global quantities, it is useful to define the mean 
profile, (u(z)). Let us consider for instance the geometry depicted in Fig. (JJJ, where the direction 
perpendicular to the walls is denoted by z. We define an homogeneous mean profile as: 



(11) 



G 




Figure 2: Results in the plane y = L y /2 along the channel measured in the transversal strip 
configuration (see left panel of Fig. QJ. The left panel shows the velocity profile. Notice that the 
pure inlet Poiseuille flow becomes an almost perfect shear free profile in the region with si = 1. 
The right panel is meant to highlight the local differences between the pure Poiseuille flow and the 
measured profiles, showing the result for (u x (r) — uP ols (r)). Notice the recirculation area, entering 
deep in the channel bulk, produced by the alternating slip and no-slip boundary conditions. 



(u(z)) = ~ J u(r) dxdy (14) 

where (■ • ■) stands for averaging over a plane parallel to the boundary surface, S. Even though 
the local velocity does not reproduce a Poiseuille profile, it can be shown from that in the 
case of periodic boundary conditions between inlet and outlet flows, the mean homogeneous profile 
(|14|l cannot develop non linear stresses, namely: 

(u(z)) = u pois (z) + u sHp (15) 

with the notable fact that a slip velocity may appear at the boundary. In the above definition, 
(|15fl u pois (z) stands for the Poiseuille parabolic profile with zero velocity at the boundary. A first 
set of qualitative results are plotted in Fig. J5J), where the local velocity profiles and the difference 
between the observed velocities and the standard no-slip Poiseuille flow are shown. 

From the expression l|15fl . one may define a macroscopic, global slip length, as the distance 
away from the wall at which the linearly extrapolated slip profile l|15|) vanishes: 

Ls = Ifl USUP < M ( 16 ) 

where \d z u po i S (z w )\ is the Poiseuille stress evaluated at the wall. Similarly, one may define the 
mass flow rate gain G as 

0.±-(l + %) (17) 



being <i> s = J u x (r)dydz the real mass flow rate and <£> p the Poiseuille mass flow rate for our 
configuration: 

/ dP\ L 3 Ly 

*»=(-^J-^f (18) 

with [i the dynamic viscosity of the fluid. In terms of these quantities one can define the 
pressure drop reduction, 
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n= AP " '^ p Ap (19) 

no— slip 

which is defined as the gain with respect to the pressure drop corresponding to a non slip chan- 
nel with the same overall throughput, <I> S . The pressu re drop reductio n, II, is usually interpreted 
as an effective drag reduction induced by the slippage. lOu et al\ l)2004|) . 

4 Numerical results 

Next, we present the numerical results obtained from the H-LBE model by changing the spatial 
distribution and intensity of the slip function at the boundaries. We shall also address dependencies 
of the slip flow on the Knudsen and Reynolds numbers. 

We begin by investigating the dependency of the macroscopic slip length, L s , and the average 
mass flow rate through the channel, on the total amount of slip material deposited on the surface. 
The natural control parameter to investigate this issue is the average of the slip function on the 
boundary wall: 



(s(r w )) = ~ Js(r w )dS (20) 



that is best interpreted as the renormalized effect of the total mass of hydrophobic material de- 
posited on the surface, at the (unresolved) microscopic level. 

Second, we also present results as a function of the non-homogeneity of the hydrophobic pat- 
tern. This non-homogeneity, or roughness, can be taken as the spatial variance of the slip function: 

A 2 = ((s(r w ) - s av f) = i J (s(r w ) - (s)) 2 dS. (21) 

In order to quantify the gain or the loss in the slip flow with respect to the homogeneous situ- 
ation, we shall focus our attention mainly on the simplest non-trivial inhomogeneous boundary 
configurations sketched in Fig. (JIJ. 

This corresponds to a periodic array of two strips. In the first strip, of length H, the slip 
coefficient is chosen as s{r w ) = S\. In the second strip (of length L — H), we impose s(r w ) = sq. 
We will distinguish the two cases when the strips are oriented longitudinally or transversally to 
the mean flow. In these configurations, the total mass s av is given by: 

Sav = £.31 + (1 - £) s 

and the degree of non-homogeneity, or roughness, by: 

A 2 = £(l-0(si-so) 2 . 

By choosing (without loss of generality) s\ > so, in this configuration the quantity £ = H/L 
is a natural measure of the localization of the slip effect. This geometry allows us to compare our 
results with some analytical, numerical and experimental results for the small Knudsen regime and 
also to extend the study to the transition regime. In sub-section (|4.3|l we shall also present results 
with slightly more complex boundary conditions, namely for the case of a bi-periodic pattern of 
alternating slip and no-slip boundary conditions. 

4.1 Exact Results and Knudsen effects 

As a validation test, we first check whether our model can reproduce some of the existing results 
concerning the slip properties of hydrodynamic systems with boundaries made up of alternating 
strips of zero-slip and infinite-slip lengths. 



iPhilid dl972^ analyzed this situation using the Navier-Stokes equations for the case of a 
cylinder with boundaries made up of alternating longitudinal strips of perfect-slip and no-slip. He 
obtained the following exact result: 

e" 9 ^ = -log(l/ COS «/2)) (22) 

Ly 7T 

where £ is the fraction of the plate where the slip length is infinite and where we have defined 
giong ag j-jjg norma ii zec j (to the pattern dimension) macroscopic slip length. 

Notice that the r.h.s. of l|22|l is independent of the radius of the cylinder, and therefore Philip's 
result is direct ly applicable to our ge ometry of Fig. Q, in the limit of small Knudsen numbers. 
More recently. lLauga fc Stonel I 20031) . analyzed the same situation with the only variant of using 
transversal rather than longitudinal strips. In the limit of a cylinder with infinite radius (plane 
wall boundaries), their result for the normalized slip length can be written as: 

e a " s ^ = -log(l/ COS «/2)). (23) 

L x 7T 

In our language, local infinite (zero) slip lengths can be obtained by choosing s\ = 1 (so = 0). 
A consistency check for our mesoscopic H-LBE model is to reproduce the hydrodynamic limits 
studied in the aforementioned papers, in the limit of small Knudsen numbers and large channel 
aspect-ratio, L z /L x . 

To this purpose, we performed a direct numerical simulation of the H-LBE model for a channel 
with square cross-section, L x = L y , and different heights, L z . For small and fixed Knudsen number, 
by increasing the aspect ratio L z /L x at fixed channel length, L x , the previous hydrodynamic 
limits are attained and the normalized slip lengths i t 8 rana ,rf n3 are independent on L z . In Fig. J3J 
we present the results obtained for both longitudinal and transversal strips compared with the 
analytical predictions (|22I23(I for a given channel aspect-ratio. 

The result (see Fig. Q) shows that the analytical hydrodynamic results are well reproduced 
by our mesosco pi c model. Moreover, w e can go beyond the hydrodynamical limit studied by 
IPhilid lll972albh : lLauga fc Stonel (2003), and investigate the effect of larger Knudsen numbers 
on these configurati ons, both in t h e nea r-hydrodynamic and in the transition regimes observed 
in the experiments, iMaurer et al\ l|2004|) . The result (see Fig. ©) is that an increase of the 
Knudsen number leads to an increase of the slip length, without preserving the ratio between 
giong an j gtrans ^ gee mse t f pjg (j*|Jj) These results can be explained by observing that upon 
increasing the Knudsen number, even the 'non-conductive' strips which had zero-slip length in 
the hydrod ynamic regime, acquire a non-zero slip due to effects of order Kn 2 in the boundary 
conditions. Isbragaglia fc Succil l)2004() . As a result, the local slip length (no longer equal to zero) 
is incremented, thereby yielding a net gain in the overall slippage of the flow. Let us notice that at 
still relatively small Knudsen numbers, Kn = 0.05, a fairly substantial increase of the slip length is 
observed, which may reach 60 — 80% of the typical pattern dimension for a percentage of slipping 
surface £ ~ 0.8. 

Another interesting question concerns the dependency of the local velocity profile on the local 
slip properties with changing Reynolds and Knudsen numbers. We choose a transversal periodic 
array of strips with H = L x /2 and sq = 0, s\ = 1 and look at the profiles in the middle of the 
region with si — 1 and in the middle of the region with sq = 0. The DNS results (Fig. Q) 
clearly indicate a dependency on the Knudsen and Reynolds numbers only in the slip region. 
This is readily understood by observing that the Reynolds number is given by Re = Ma/Kn, so 
that, by fixing the Mach number and varying the Reynolds number, we change also the Knudsen 
number, thus affecting the local slip properties of the flow. The most interesting result here is the 
inversion of concavity for the local profile nearby the wall in the slip region: a clear indication of 
the departure from the parabolic shape of the Poiseuille flow. 

Next, we check our method against some experimental results and MD simulations. For ex- 
ample, in Fig. (JjjJ, we show the dependency of the transversal normalized slip length, £ t J' ans , as a 
function of the inverse of the Knudsen number, i.e. as a function of the mean channel pressure, 
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Figure 3: Normalized slip length for transversal and longitudinal strips with s\ = 1, sq = 0. We 
plot the normalized slip length as a function of the slip percentage £. The system's dimensions 
are those of Fig. (QJ. A first set of LBE simulation is carried out at small Knudsen, Kn = 1.10 -3 
for transver sal (□) and lo ngitudinal strips ( o). These results are compared with the analytical 
estimates of lPhilid l)l972a|) (dashed line) and lLauga k. Stone] l|2003l) (continuous line). Notice the 
perfect agreement with the analytical results in the hydrodynamic limit. Another set of simulations 
is carried out with much larger Knudsen, Kn — 5.10~ 2 to highlight the effect of rarefaction on the 
system for both traversal (x) and longitudinal (+) strips. In the inset we show the ratio between 
the slip lengths for parallel and longitudinal strips for Kn = 1.1CP 3 (o) and Kn = 5.10 -2 (□). 
Here we notice how by increasing the Knudsen number the orientation of the strip region with 
respect to the mean flow becomes less important. 
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Figure 4: Local velocity profile in the middle of a slip strip (the one with si = 1) for transversal 
strips in the geometry depicted in Fig. JJ). We plot the velocity in the stream- wise direction as 
a function of the height z for two different Reynolds numbers Re ~ 4.5 (□), and Re ~ 9.5 (o). 
The Re numbers are estimated as the ratio between the center channel velocity of the integral 
profile and the sound speed velocity c s . Both velocity fields are normalized with the center channel 
velocity. The Knudsen numbers are Kn — 0.01,0.005 respectively. Inset: the same but in the 
middle of a no-slip strip (s = 0) 
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Figure 5: Left panel: Normalized transversal slip length, £* rans = as a function of the average 
pressure in the system (inverse of the Knudsen number) and for different values of the localization 
parameter: £ = 0.58 (□), £ = 0.65 (o), £ = 0.63 (A). The values of so,si are kept fixed to 
so = and s\ = 1. This behavior is qualitatively similar to what observed in MD simulations of 
microchannels with grooves, where the degre e of slippage localization is governed by the width of 
the grooves, see ICottin-Bizonne et all (120041) . In the inset we show the same trends but for the 
mass flow rate gain G (see eq. I17fl . Right panel: £ t s rans as a function of the local degree of slippage 
so for different values of slippage localization, £ = 0.25,0.5,0.75 (D,o,A respectively). In the inset 
we show the dependency of the slip length, £* rans ; on the microscopic slip properties, sq/(1 — sq), 
for the same values of £. 



for differ ent values of the localization parameter £. This is a direct comparison with the results in 
Fig. 6 of ICottin-Bizonne et all l|2004j) where the evolution of the slip length as a function of the 
Pressure in MD simulations of a channel with grooves of different width is shown. Also in that 
case, the slip length increases by either decreasing the pressure (increasing Knudsen) or increasing 
the groove width (increasing the region with infinite slip). The two behaviors are qualitatively 
similar, with a less pronounced slip length for our case also due to the fact that we show the 
case of transversal strips while in Fig. 6 of ICottin-Bizonne et all l|2004h only the case of longi- 
tudinal grooves are presented. In the right panel of fig. (JHJ we plot £ t s rans at varying the level 
of slippage, so, of one of the two strips (the other being kept fixed to s\ — 1). This is meant to 
investigate the sensitivity of the macroscopic observable to the microscopic details. As one can 
see, the change is never dramatic, at least for this configuration. In the inset of the same figure 
one notice a linear dependency between £ t s rans j and the local slip properties, so/(l — So). For 
local slip properties we mean the local slip length as defined from the local boundary condition, 
u\\(r w ) oc ^( r "') |d n ui|(r m )|. The same linear trend is observed in fig. 12 of lCottin-Bizonne et all 
(2004) using a hydrodynamic model with suitable boundary conditions. 

I n Fig. we pre sent the same kind of plot shown in the experimental investigation (see Fig. 
15 of lOu etcdl l|2004|) '). Here, we plot the pressure drop reduction, II, in the microchannel as a 
function of the percentage, £, of the free slip area on the surface (super-hydrophobic material). We 
note a remarkable agreement with the experimental results over a wide range of £, i.e. the ratio 
between the regions with super-hydrophobic and normal material on the wall. The geometry and 
Knudsen number are the same of the experiment. The only free parameters are the values of So and 
Sj assigned to the two different strips. Here, we have fixed si = 1 in the super-hydrophobic area, 
and we have varied so S [0.4, .65] for the normal material. Notice that the LBE results exhibit the 
same trend of the experiments as a function of £ and they are even in good quantitative agreement 
for £ ~ 0.6. Overall, there is a small dependency of n on the unknown value of so, at least in the 
range considered, as already shown by the data presented in Fig. J^J. Once the correct values of 
si and so able to reproduce the experimental results are identified, one may easily use the present 
LBE method to predict and extend the outcomes of other experiments with different geometries 
and/or distributions of the same hydrophobic material on the surface. 
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Figure 6: Results for the pressure drop reduction IT Ijl9|l . as a function of the percentage of the 
free slip area on the surface. The geometry is the same investigated experimentally in Fig. 15 of 
lOu et al\ l)2004[) where a micro-channel with only one surface engraved with longitudinal strips 
with super- hydrophobic material is studied (left panel). Here we present the results from the 
experiments (•), superposed with our numerical simulations (x) obtained at Kn = 5.10 . The 
numerical mesh is such to mimick the same channel height of 127/^m used in the experiment. The 
error bars in the LBE numerics represent the maximum variation obtained by fixing s\ = 1, in 
the shear-free area, and changing sq G [0.4, 0.65] in the normal surface. 



4.2 Effects of roughness and of localization 

As a next step, we address of the effect of the roughness, the total mass of the hydrophobic 
material, and the localization over the surface. 

To this purpose, we choose a transversal configuration where H = L x /2 (fixed localization, 
£ = 0.5) and 

so = s av + A, si = s av - A, 

thus yielding (s) = s av for any degree of the roughness, A. 

We then look, for a given s av , at how the slip properties of the system respond to changes in the 
excursion, A, at fixed localization £. The results for the normalized transversal slip length, /' mns , 
and the mass flow rate gain, are presented in (Fig. (JJJ)). Both the slip length and the mass flow 
rate increase by increasing A. Notice that one can easily reach slip lengths which are of the order 
of 10% of the channel pattern dimension by increasing the roughness at fixed total mass of slip 
material deposited on the surface. Similarly, the mass flow rate gain, G, is increased of the order 
of 20 — 30% with respect to the Poiseuille flow. In other words, the best throughput is obtained 
by increasing the inhomogeneity of the slippage material deposited on the surface (experimentally 
this means to keep the region covered by hydrophobic molecules as segregated as possible from 
the region covered with hydrophilic molecules). 

Another possible way to compute the slippage effects is to analyze the slippage at fixed total 
mass of hydrophobic material, varying both the localization and the roughness. To this purpose, we 
choose again a transversal configuration where we fix the total mass s av and we choose s% — s a v/£, 
sq = 0, £ = H/L x being the degree of localization associated with a given roughness. The result 
(Fig. (2)) is that the slippage is greater as the the degree of localization, and -consequently- of 
roughness, is increased. 

4.3 Mean field approach and beyond 

We observe that a mean field approach is able to reproduce the qualitative trends observed so far. 

In the boundary condition (|12|l there is a coupling between the local velocity field and the 
stresses at the wall. Obviously, the averaged slip length depends both on s(r w ) and on the 
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Figure 7: Normalized transversal slip length, £ t s rans = L s /L x , for boundary conditions as in Fig. 
JTJ with si — s av + A in a fraction £ = H/L x = 1/2 and so — s av — A in the other. The Knudsen 
number is Kn = 0.001. We plot the normalized slip length as a function of the roughness parameter 
A and the following values of s av are considered: s au =0.4 (□), s av — 0.5 (o), s av — 0.6 (A). Inset: 
Same trends of the main body but for the mass flow rate gain G. 
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Figure 8: Normalized transversal slip length, £* ra " s = L s /L x , for boundary conditions of transver- 
sal strips as in Fig. with si = s a v/ti in a fraction ^ and so = in the other. The Knudsen 
number is Kn = 0.001. We plot the normalized slip length as a function of the localization pa- 
rameter £ for the following values of s a v- s av — 0.25 (□), s av — 0.5 (o), s av = 0.75 (A). Inset: 
Same trends of the main body but for the mass flow rate gain G. 
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stresses. 

In order to highlight the effect of the slip function on the mean quantities, as a first approxi- 
mation, we can leave the wall stress fixed at its Poiseuille value, and work only on the properties 
of s(r w ), namely: 

(usU P ) ~ (y^>. (24) 
For the configuration analyzed so far, without loss of generality we define: 

s- si - s s + si 

O = ^ Si = 

2 + 2 

and write the averaged slip properties (jzr^) as a function of 8 and s + : 

/ s s si s+-5 s++5 

-)=Po-, 1- Pi- =Poz — F+Pii ? I 25 ) 



1 — s 1 — so 1 — si 1 — s+ + 8 1 — s+ — 8 

where p\ = jf-, po — 1 — Pi are the percentages of the surface associated with slip and no-slip 
areas respectively. Making use of Taylor expansion up to second order in 8 we get: 

; s . s + 8 , . S 2 .„„. 

i ^r^+n ufo-lxO + Ti 26 

1 — S 1 — 3+ (1 — s+) z (1 — S+) J 

Since Pi = jf- = £, and po +Pi = lj we finally obtain: 

' S ) M S+ | 8(2£-l)(l- s+ ) + 6^ (2?) 



x l-s' l-s+ (l-s+) 3 

First, in our case of a fixed localization, by setting £ = 1/2 we have s+ = s „, 8 = A and we 
obtain: 

/ 5 \ ^ Sat) _|_ ^ ^8) 
1 s 1 s av (1 s a ^) 3 

that results in a greater slippage when the roughness A is increased. 

Second, if we choose si = and so = 0, as for the case with fixed total mass, we obtain 
8 = 4gp and s + = 4gp ■ This results in 

/ s ^ s a . u Sa „(4g 2 -2Q(2^- Sa „) + 2^ 

(2£-a „) + (2e-s OT ) 3 1 j 

that, as a function of the localization ^ yields a qualitative agreement with our analysis, supporting 
the idea that the effect of slippage is greater when slip properties are localized. 

It should be appreciated that the mean field approach discussed above is not exhaustive. In fact, 
we can design an experiment with the boundary configuration sketched in Fig. , and investigate 
the total slippage as a function of the distance, d, between the strips. For this geometry, the mean 
field approach presented before would yield the same results irrespectively of d. 

On the other hand, we expect non-linear effects to be present when the strips get close enough, 
due to the correlation between s(r w ) and the stress at the boundary, d n u(r w ). Indeed, as one can 
see in Fig. Jj]J, the slippage is increased when the two strips get closer to each other. 

This effect, even if only of the order of 10% in the mass flow rate with respect to the con- 
figuration for d > 1, cannot be captured by the previous mean field argument. Let us notice 
that a similar sensitivity to the geometrical pattern of t he slip and no-slip areas has been recently 
reported in the experimental investigation of lOu et~al\ l)2004|) . where it is found that for the same 
microchannel geometries and shear-free area ratios, microridges aligned in the flow direction con- 
sis tently outperform re gular arrays of microposts. Similar considerations have been also presented 
bv IVinogranoval 1 1999). 
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Figure 9: Left: The configuration with a transversal bi-strip structure at the walls. The total 
boundary lengths are L x , L y (stream, span). The slip coefficient is chosen as si = 1 in two strip of 
length H and s = in the others. The distance between the two strips is d. Periodic boundary 
conditions are always assumed in the span- wise and stream- wise direction. Right: results for the 
slip length and the mass flow rate gain (inset) as a function of the distance d (in lattice units) 
between the two free-shear (s = s\) strips of width H = 20 (in lattice units). All the other 
parameters, L x = L y = 64, L z = 84, Kn = 1.10~ 3 are kept fixed. 

5 Conclusions 

We have presented a mesoscopic model of the fluid-wall interactions which proves capable of 
reproducing some properties of flows in microchannels. We have defined a suitable implementation 
of the boundary conditions in a lattice version of the Boltzmann equation describing a single-phase 
fluid in a microchannel with heterogeneous slippage properties on the surface. In particular, we 
have shown that it is sufficient to introduce a slip function, < s(r w ) < 1, defining the local degree 
of slip of mesoscopic molecules at the surface, to reproduce qualitatively and, in some cases, even 
quantitatively, the trends observed either in MD simulations or is some experiments. The function 
s(r w ) plays the role of a renormalizing factor, which incorporates microscopic effects within the 
mesoscopic description. 

We have analyzed slip properties in terms of slip length, L s , slip velocity, V s , pressure drop 
gain, II and mass flow rate $ s , as a function of the degree of slippage, and its spatial localization. 
The latter parameter mimicking the degree of roughness of the ultrahydrophobic material in real 
experiments. 

With a proper choice of the slip function s(r w ) in longitudinal and transversal configurations, 
we have reproduced previous analytical results concerning pressure-driven hydrodynamic flows 
with boundarie s made up of alternatin g strips of zero-slip and infinite-slip (free-shear) lengths, 
IPhilid lll972albh : lLauga fc Stone! lE)0.T). We have also discussed the increment of the slip length 
in the transition regime, i.e. where the Maxwell-like slip boundary conditions Q are supposed to 
be replaced by second-order ones (J2J • 

The local velocity profile has also been studied with changing Reynolds and Knudscn numbers 
and the local slip properties on the surface. 

The method introduced is able to describe slip lengths of the order of the total height of the 
channel (of the order of tens of /Ltm), or fractions thereof. This is accompanied by an important 
increase in the mass flow rate, or equivalently, in the pressure drop gain. Whenever possible, we 
have compared the results based on the Heterogeneous LBE with MD simulations and with some 
recent experiments. 

In particular, we have shown that the H-LBE approach is able to reproduce the increase of 
the slip length as a functi on of the inverse of the mea n pressure in the channel, as observed in 
recent MD simulations bv ICottin-Bizonne et all l)2004h . Concerning the same MD simulations, 
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Figure 10: 2d and 3d lattice discretization of the velocities in the LBE schemes used in this work. 
The velocities entering in the north wall in the 2d scheme (left) are, f5,f2,fe- The outgoing 
velocities are: fy, fi, /g. At the south wall the roles are exchanged. 

we have found a similar linear dependency of the macroscopic slip lengths, L s as a function of 
the microscopic slip properties at the surface. As to the experiments, we have shown that the H- 
LB E app r oach i s able to achieve quantitative agreement with the experimental study presented in 
lOu et an l|2004|) . concerning the slip properties as a function of the relative importance of regions 
with high-slip and low-slip at the surface. The natural application of our numerical tool consists 
in tuning the free parameters so and s\ in order to reproduce experimental results in controlled 
geometries. 

Then, one may use the LBE scheme with the given sq and s\ values, to explore flows in different 
geometries and/or with different patterns of the same slip and no-slip materials. 

The method is a natural candidate to study flow properties in more complex geometries, of 
direct interest for applications. Transport and mixing of active or passive quantities (macro- 
molecules, polymers etc..) can also be addressed. 

By definition, the present H-LBE description is limited to a phenomenological interpretation 
of the slip function. A natural development of this approach, is to implement a multi-phase 
Boltzmann description, able to attack the wall-fluid interactions and fluid-fluid interactions at a 
more microscopic level. 

This route should open the possibility to discuss the formation of a gas phase close to the wall, 
induced by the microscopic details of the fluid- wall physics. Results along this direction, will make 
the object of a forthcoming publication. iBenzi et al\ l|2005f) . 



6 Appendix A 



The Lattice Boltzmann Equation (LBE) for a Pressure-driven channel flow is a streaming and col- 
lide equation involving the particle distribution function // (x 1 1) of finding a particle with velocity 
q (discrete velocity phase space) in x at time t. The equation is written in the following form: 

f l (x + c l S t ,t + S t ) - /,(*,*) = ~ (/,(£, i) - rt eq) (p,itj) + 5 ^F 9l (30) 
with r the relaxation time and gi the forcing term projection with the property 

5> = 5><3 = L (31) 

/ i 

For the case of two dimensional grid (2DQ9) depicted in Fig. (|10|l . for example, the gi's can be 
taken with the following properties: 

9i = ~93 55 = 98 = -.96 = -97 (32) 
leaving only one unknown parameter, say 55. Discrete space and time increments are S Xl 5t, with 



St 



the intrinsic lattice velocity. The equilibrium distribution f^ eq ^ (p, u) is given by: 



fi (p>rf) = w ip 



(q ■ u) 1 (cj ■ u) 2 _ 1 
~~ 2~ 3 2Z 



(33) 
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being c 2 s — -|c 2 the sound speed velocity. Concerning the 2DQ9 model here used for the technical 
details, the velocity phase space is identified by the following discrete set of velocities: 

c = (0,0)c 

ci,c 2 ,c 3 ,c 4 = (l,0)c,(0,l)c,(-l,0)c,(0,-l)c (34) 
c 5 ,c 6 ,c 7 ,c 8 = (1, l)c, (-1, l)c, (-1, -l)c, (1, -l)c 

and the equilibrium weights are w = 4/9, u>; = 1/9 for I — 1, ...,4, = 1/36 for Z = 5, ...,8. 

As far as the 3d model we use in the numerical analysis (3DQ19), it is a 19 velocity model whose 

velocity phase space is identified by: 

c = (0,0, 0)c 

ci,2,C3.4,c 5 , 6 = (±l,0,0)c,(0,±l,0)c,(0,0,±l)c (35) 

C 7 ,...,10, Cn,...,i4, Ci5,...,18 = (±1, ±1, 0)C, (±1, 0, ±1)C, (0, ±1, ±1)C 

and equilibrium weights u>q = 1/3, u>i = 1/18 for Z = 1, ...6, wi = 1/36, I = 8, 19. 

Our hydrodynamic variables such as density p and momentum pu are moments of the distribution 

function f t = f t (x, t): 

i i 

and in order to derive Hydrodynamic equations from H3U[) . we must consider the following expan- 
sions: 

nl 



f l {x + c l 8 u t + 8 t ) = Y^—D\ l f l {x 1 t) (37) 

OO 

fi = E e "// n) ( 38 ) 



n=0 

OO 



9 * = E ( 39 ) 



n=0 

where e — 8 t and D t = (9 t + ci ■ V). order by order in e and l|3U|l imply: We can use the expansions 
(|37|I . (|38() . H39|I in (|30() and by equating order- by-order in e self-consistent constraints on fj are 
obtained. 

Up to the first order in e with dt = dt + edt ± we obtain the following equations: 



d tP + V{pu) = 

d t u + (u • V)u = -~VP + ±V • (i/pVu) 



(40) 



with = ^ T 3 2 ^ ^ and where VP contains both the imposed mean pressure drop, F, and the fluid 
pressure fluctuations. 



7 Appendix B 

Let us now go back to eq. (|30|l . and derive explicitly the non- homogeneous boundary conditions 
used in the text l|12l) . in the limit 8 X = St — > 0, with c = 4* — * 1. We specialize to the steady-state 
boundary condition at the north-wall (z = L z ) for the 2d lattice (2DQ9): 

h{^w) = (1 - s(r w ))f 5 (r w ) + s(r w )/ 6 (r w ) 
fa(r w ) = f2(r w ) 

fa(r w ) = (1 - s(r w ))/ 6 (r w ) + s(r w )/ 5 (r TO ). 
Assuming a constant density profile p = 1 in the fluid, by definition we have for r = r w : 

u \\( r w) = fi(r w ) " fa(r w ) + /s(r w ) - /e(r^) + / 8 (r m ) - h{v w ). (41) 
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In the limit of small Mach numbers, disregarding all 0(u 2 ) terms in the equilibrium distribution 
and using the steady state, dtf = 0, expansion: 

oo 

/,(r) = ^(-l)"(r(c ; V)n// e9) (p,u) +rFt] (42) 

n=0 

we finally obtain the estimate for the slip velocity Uii: 
u\\(r w ) = 2Frg 1 + %u\\(t w ) + ~c 2 T 2 d x u\\(r w )+ 

+2s[2Ft 95 + _ fd x u ± (v w ) - fd y u n (r w ) + ^f(d 2 x + 8$)u\\(t w ) + ^fd x d y u {l (r w )] + 0(t 3 ). 

(43) 

By noticing that the external forcing, F, is of the order of magnitude of the second-order stress, 
\dyU\\(Y w )\v, and v = c 2 t, the first order in Kn of 143|1 reads: 

u N (r„) = Kn (£) 1 £M^ ]9 wU ||(r w )| (44) 

where we have used L z d y — d n , r = L ^ n ; ^ui^r^) = — |9 n ui|(r w )|, which is the expression 
used in the text. 

The second-order calculation in Kn is particularly simple if we specialize to an homogeneous 
case (<9 X (») = 0). After some calculations, we obtain: 

^=Kn(^^ s \d n ^\+Kn 2 (^j (1 - 4 ffB )|^U|,| (45) 

which is the secon d order, in Kn, boundary conditions, with unknown parameters s and 55 (0 < 
.95 < 1/4 ) used bv lSbraeaglia fc Succil l|2004 - 
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